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Abstract 



The quantum density of states of the Henon-Heiles potential displays a 
pronounced beating pattern. This has been explained by the interference 
of three isolated classical periodic orbits with nearby actions and periods. 
A singular magnetic flux line, passing through the origin, drastically alters 
the beats even though the classical Lagrangean equations of motion remain 
unchanged. Some of the changes can be easily understood in terms of the 
Aharonov-Bohm effect. However, we find that the standard periodic orbit 
theory does not reproduce the diffraction-like quantum effects on those clas- 
sical orbits which intersect the singular flux line, and argue that corrections 
of relative order Ii are necessary to describe these effects. We also discuss the 
changes in the distribution of nearest-neighbour spacings in the eigenvalue 

spectrum, brought about by the flux line. 

PACS numbers: 05.45. +b, 02.50.+s,03.65.-w,03.65.Ge 

Typeset using REVTgX 

1 



I. INTRODUCTION 



The connection between the quantum behaviour of a particle and its classical motion 
in a given potential is of continuing interest. Particularly revealing is the link between the 
oscillating part of the quantum density of states (obtained by subtracting out its smooth, 
Thomas-Fermi like component) and the classical periodic orbits [1,2]. Orbits with short 
periods govern the large-scale energy dependence of the quantum density, while orbits with 
long periods determine its fine structure. It is known that even when the classical motion of 
a particle is globally chaotic, isolated periodic orbits may continue to exist. For that case, 
Gutzwiller [1] has derived the so-called "trace formula" for the oscillating part of the density 
of states: 
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Sg(E) = J2Y,Axk(E) cos kl-Sx(E) 
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Here A goes over all primitive periodic orbits (supposed to be isolated here); k counts the 
number of revolutions around each primitive orbit, yielding a series of harmonics; S\(E) = 
§p\ ■ dq\ is the classical action integral along the primitive orbit A; and the Maslov index 
(T\ is a phase depending on its topology. The amplitudes A\k are given by 
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in terms of the period T\ = dS\/dE and the stability matrix M\ of the corresponding orbit 
(7 being the 2x2 unit matrix). 

If the fine structure in the quantum density of states is erased by an appropriate smooth- 
ing [2-4], then the remaining oscillations may be attributed to the lowest harmonics of the 
periodic orbits with the shortest periods. The objective in this paper is to study this link 
between the orbits of shortest period and the "coarse-grained" quantum density of states, 
particularly in the presence of a magnetic flux line. 

The potential of Henon and Heiles (HH) [5] is very appropriate for such a study. A 
classical particle can escape from this potential over three barriers if its energy is larger 



than the threshold energy. Below threshold the motion is confined; it becomes less regular 
with increasing energy and is fully chaotic at and above threshold [5]. If the particle energy 
is less than some 80 percent of the barrier height, then it is established [6,7] that there are 
only three types of distinct isolated periodic orbits of nearly equal periods and actions, all 
others having periods about twice as large or more. As we have shown in our earlier paper 
[7], the interference of these three periodic orbits can explain the beats in the quantum 
density of states. 

In the present work, a magnetic flux line through the origin is added to the HH potential. 
Because of the beats in the quantum density of states, the HH potential provides a new 
testing ground to study the effects of a flux line and their connection to classical orbits. We 
find that the beat structure in the quantum density is altered drastically in the presence 
of the flux line. Moreover, the vortex potential which is singular at the origin introduces 
an additional zero in the s-state components of the quantum wavefunctions. The Fourier 
analysis of the quantum level density reveals, indeed, a new peak of approximately half 
the original period. In view of the fact that the classical Lagrangean equations of motion 
remain unaltered by the introduction of the flux line, the question arises if the changes in 
the quantum level density may be explained within the standard periodic orbit theory by 
simply adding the Aharonov-Bohm phase to the actions of those orbits that have a non-zero 
winding number about the origin. The answer is found to be negative. 

The nearest-neighbour spacing (NNS) distributions of the quantum spectra show very 
interesting changes with the introduction of the flux line. The HH potential has a three-fold 
symmetry, and the quantum states can be classified in three distinct classes. Two of these 
classes have identical spectra in the absence of the flux line. This degeneracy is broken 
by the latter, and nontrivial changes in the NNS distributions take place. In particular, 
with the introduction of the flux line, the distribution for one symmetry class changes from 
Wigner to Poisson type - a purely quantum-mechanical effect. 

Our paper is organized as follows. In section II, the quantum beats of the HH spectrum 
(without magnetic flux line) and their semiclassical interpretation are described. The content 

3 



of this section is close to that of Ref. [7], except for a more careful calculation of the 
Maslov indices following the methods of Creagh, Robbins and Littlejohn [8]. Moreover, 
the numerically determined actions of the three leading orbits have been parametrized more 
accurately. We also re-analyze their degeneracy factors due to the symmetries of the HH 
potential, and find that Eq. (1) reproduces the quantum level density much better if one of 
the orbits is given an extra factor of three. 

In section III, a singular flux line through the origin is introduced. We first describe the 
quantum results in detail and then try to interpret them in terms of the standard periodic 
orbit theory. A simple modification of the classical calculation, just shifting the action of 
the loop orbit that encloses the flux line by the Aharonov-Bohm phase, is not sufficient to 
fully explain the quantum results. In section IV, the NNS distributions are displayed and 
discussed in detail. Section V gives a short summary and the Appendix contains an account 
of the Maslov index calculation without the flux line. 

II. BEATS IN THE HENON-HEILES DENSITY OF STATES 

The HH potential [5] for the planar motion of a particle is given by 



Following the usual convention, we put h=m=u=l, and a is a dimensionless parameter 
whose strength determines the anharmonicity. The equipotentials of Vhh are shown in 
Fig. 1; they have a three-fold symmetry. This is obvious by writing Vhh m polar coordinates 
(r, 6) in which it has the form 



Along the three symmetry axes 6 = 7r/2, (tt/2 + 27r/3), and (7r/2 + 47r/3), we have sin(3#) = 
— 1 so that the barrier height is minimum; the threshold energy there is E* = 1 /6a 2 (in 
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We have obtained the quantum spectrum of the Hamiltonian 

H = \(PI+P 2 y ) + V HH (5) 

by its diagonalization in a large harmonic oscillator basis \nl) with radial quantum number 
n and angular momentum / (and oscillator constant u = 1). Because of the noncentral 
nature of Vhh, I 1S n °t a good quantum number and states with A/ = ±3 are mixed. This 
mixing takes place within three distinct and disconnected sets or symmetry classes of the 
basis states: 

{/}: / e {..._6, -3, 0, 3, 6...}, 
{//}: Ze{...-5, -2, 1, 4, 7...}, 
{///}: Ze{...-4, -1, 2, 5, 8...}. (6) 

which together cover the full Hilbert space. Note that under time reversal, /—»—/, and 
therefore set {/} maps onto itself, while {//} maps onto and vice versa. In the 

absence of an external magnetic field, the Hamiltonian H of Eq. (5) is invariant under time 
reversal, and thus the sets {//} and yield identical eigenspectra. Thus H may be 

diagonalized separately in each basis set. Another interesting point to note is that only 
the states belonging to the set {/} are invariant under the rotations 6 — > (6 + 2tt/3) and 
6 — > (6 + 47r/3). Therefore only the states of set {/} have the same three-fold symmetry 
as the Hamiltonian, but not the states belonging to the sets {//} and {///}. This will be 
important for the analysis of the NNS distributions in the different spectra, which will be 
done separately for each symmetry class. The quantum spectrum obtained by diagonalizing 
H is discrete even for E > E* (the barrier height) due to the truncation of the basis, and 
is only reliable for E < E*. The basis size used in the present calculations has an energy 
cut-off E cut = 130 ifiu). It should be noted that in principle, even the "bound" states 
below the threshold have non-zero widths and the particle has a finite probability to tunnel 
through the barrier. The time scale relevant to our analysis is certainly much smaller than 
the tunneling times, so that the neglect of these widths should not affect our conclusions. 
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As in our earlier paper, a Strutinsky smoothing [3] is applied to obtain the overall smooth 
part of the density of states. Let us write the exact quantum-mechanical level density as 

g(E) = J2HE-E t ), (7) 

i 

where the sum over {i} includes all three sets of eigenstates. To obtain the Strutinsky 
smoothed density g (E) (which is equivalent to the extended Thomas-Fermi value [9,10]), 
each delta function in Eq. (7) is replaced by a Gaussian of half- width 7 > 1, multiplied by a 
so-called "curvature-correction" function [3] that may be written as a generalized Laguerre 
polynomial: 



g(E) = ^-=Y^ ex V 
7VT i 



(E - Ei) 
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It can be shown [3,9] that the inclusion of this correction polynomial preserves the 2m lowest 
terms of a local Taylor expansion of the average level density. Here we use a sixth-order 
polynomial (m=3). For 7 ~ 1.4, all oscillations in the level density are wiped out and the 
so-called "plateau condition" (local independence of g (E) of 7) is well fulfilled. Next, a 
Gaussian smoothing (without curvature-correction) is performed with ^ osc < 0.6, in order 
to obtain a "coarse-grained" quantum density of states which we denote by g losc (E)- Note 
that this is equivalent to suppressing all higher harmonics (k > 1) in the Gutzwiller trace 
formula (1) as well as the contributions of orbits with long time periods. Our coarse-grained 
oscillating density of states is thus defined as 

6g(E) = [g losc (E)-g(E)]. (9) 

The Sg(E) obtained by this procedure is shown in Figs. 2 and 3 (upper row) for three 
values of the anharmonicity parameter a and is essentially the same as that given in Ref . [7] . 
Note that the beating pattern for different values of a is approximately invariant when 
plotted as a function of aE. The reason for this scaling behaviour will become clear from 
the semiclassical analysis below. A discrete fast Fourier transform (FFT) of Sg(E) for 
a = 0.04 in the time domain clearly reveals three spikes around the period T ~ 1 (see Fig. 4 
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below). These may be identified with the periods of the three classical orbits A, B, and C 
that are shown in Fig. f , averaged over the energy interval (0 < E < 40) used in the Fourier 
transform. 

For obtaining the "Gutzwiller connection" to the quantum beats by means of the trace 
formula (1), it is essential to calculate the classical actions S\ of the three periodic orbits 
A = A, B, and C. This is done from the numerical solutions of the classical equations of 
motion. Note that the strength parameter a may be scaled away from the Hamiltonian H 
in Eq. (5) by defining the scaled variables £ = ax and rj = ay. The classical equations of 
motion are then universal in these coordinates with a dimensionless scaled energy e defined 
by 

e = E/E* = 6a 2 E. (10) 

The numerically calculated scaled actions s\ = 6a 2 S\ of the above three orbits can be 
parametrized, after performing a least-chi-squared fit, as 

s A (e) = 27r(e + 0.0821e 2 + 0.0585e 6 ) (0 < e < 1), 

s B (e) = 27r(e + 0.0698e 2 - 0.0046e 4 ) (0 < e < 2), 

s c (e) = 27r(e - 0.0234e 2 + O.OOlle 4 ) (0 < e < 2). (11) 

[Unlike the actions of the orbits B and C which were computed numerically, the action 
of orbit A can be written analytically in terms of an elliptic integral.] Rewriting the true 
actions S\(E) as functions of the true energy variable E, we find from the expansions (11) 
the following a dependence: 

S\(E, a) = 2tt[E + c x (aE) 2 + 0{a 6 E 4 )} (12) 

with constant numerical factors c\. The first term in the action of each orbit is a independent 
and equal to the action 2ttE of the harmonic oscillator, thus guaranteeing the correct a — > 
limit of the above fits [11]. It is the second terms, proportional to (aE) 2 , that determine 
the beat pattern for not too large values of a and E [12]. This is the reason why the energy 
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variation of the beating amplitudes of Sg(E) depends on the combination (ctE) and thus 
scales with a -1 , as can be observed in Fig. 2 when comparing the results for the different 
values of a. 

We now write down the (semi)classical expression for the oscillating part of the density of 
states, denoting it by 6g c i(E) } by just keeping the lowest harmonics (k = 1) in the Gutzwiller 
trace formula (f): 



Sg d (E) = 3A A cos S a - cr A ~ + 3A B cos S B - cr B - + 2 Ac cos S c - cr c - . (13) 



Due to the three-fold symmetry of the HH potential, there exist three distinct and energet- 
ically degenerate orbits of type A and B, whereas the orbit C maps onto itself under the 
rotations about 120 and 240 degrees. Therefore, the amplitudes of A and B have been mul- 
tiplied by a factor 3 in Eq. (13). The factor 2 multiplying the contribution of orbit C is to 
account for its two topologically different orientations, orbiting clockwise and anti-clockwise 
around the origin. 

As is evident from Fig. 3, the agreement between the quantum and semiclassical results 
is much better if the contribution of orbit C is also multiplied with a degeneracy factor 3 
(as has already been done in Ref. [7]). At present, we have no explanation for this fact, 
except for some numerical evidence from the Fourier spectrum of the quantum level density 
discussed in Sect. III. A below. If the orbit C had not the same discrete rotational symmetry 
as the HH potential, this extra factor of 3 would be justified. We find, however, no deviation 
from this symmetry within the numerical accuracy of our solutions, in agreement with the 
existing literature [6,13]. 

Besides the actions which we have already mentioned, the amplitudes A\ (2) and the 
Maslov indices a\ in (1) have been obtained numerically from the monodromy matrices 
of the three orbits (see Appendix). As already mentioned in Ref. [7], the amplitudes (2) 
diverge in the limit e — > 0, the eigenvalues of M\ becoming equal to unity for all orbits. This 
happens because in the harmonic-oscillator limit a = 0, the trace formula (1) does not apply 
since the periodic orbits of the harmonic oscillator are not isolated. We therefore regularize 
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the amplitudes, imposing their exact values 2E [see Eq. (23) below] for the 2-dimensional 
harmonic oscillator by a spline fit in the low-energy region [7,12]. In Ref. [7], this was done 
only for the orbit B and the other amplitudes were determined through the ratio of the 
unsplined Gutzwiller amplitudes. We now do the splines independently for all three orbits 
A, B, and C. 

The Maslov indices a\ were obtained in Ref. [7] simply by counting the conjugate points 
of the respective orbits. As pointed out by Creagh et al. [8], this is strictly correct only in 
one dimension, while in higher dimensions there is another contribution to the Maslov index 
arising from the saddle-point approximation in taking the trace of the energy-dependent 
Green's function. We have used the methods of Creagh et al. [8], described in the appendix, 
to calculate the Maslov indices. We find that their values for the orbits A, B, and C are 
(7^=5, <tb=4, and <7c=3, respectively. Note that these are different from Ref. [7] where their 
erroneous values and the slightly inaccurate parametrization of the actions for Sb and Sc - 
together with their different splines - appear to have compensated each other to give a fair 
agreement. 

The present calculation of 8g c i using Eq. (13) with the new Maslov indices yields a 
reasonable agreement with the quantum calculation of Eq. (9), as seen in Fig. 2 where the 
semiclassical results obtained according to Eq. (13) are shown in the lower row. The beat 
structure is reproduced with minima approximately at the correct energies. However, the 
interference is evidently not strong enough to give sufficient destruction at the beat minima. 
The disagreement is worst at the first minimum which is due the interference of orbit C with 
orbits A and B (having nearly identical periods, actions and amplitudes). 

If we increase the amplitude of orbit C by an extra factor 3 in Eq. (13), the agreement 
becomes much better and now is almost perfect. The corresponding result is shown in 
Fig. 3. Although we cannot theoretically justify this choice, it receives strong support from 
the Fourier spectra discussed in the next section (see Fig. 4). 

In passing, we remark that the a -1 scaling of the beating amplitudes in 6g(E), which 
we have explained semiclassically through Eqs. (12) and (13) above, can also be understood 
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purely quantum-mechanically. In fact, if one treats the anharmonic part of the HH potential 
(3) in second-order perturbation theory, one finds for its eigenvalues the following expression 
[12]: 



Ei(a) = E, 



(o) 

nl 



5 



+ ..., (14) 
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where E$ are the unperturbed harmonic-oscillator levels. Thus, if we take the E$ as 
the energy scale of the quantum-mechanical spectrum, we see again that the dominating a 
dependence of the approximate E{ scales with a -1 . 

III. EFFECTS OF THE MAGNETIC FLUX-LINE ON THE LEVEL DENSITY 

We now discuss the modifications brought about by the addition of a magnetic flux line 
perpendicular to the x — y plane at the origin of the HH potential. Its influence is most 
compactly expressed in the Lagrangean formalism. The Lagrangean is then given in polar 
coordinates by 

L = \(r 2 + r 2 9 2 ) - V HH (r, 9) + 86, (15) 

where 8 is the (dimensionless) strength of the flux line. The problem may also be described 
in terms of a particle moving in the HH potential in the presence of a magnetic field which 
has a delta function singularity at the origin and is zero everywhere else. The constant 
flux then is given by 8. For 8 = 1, one unit of flux quantum is enclosed by any path that 
encloses the origin. For studying the mathematical problem, we shall assume that 8 can take 
a continuous range of values between and 1, thus allowing fractional flux quanta. This is 
only allowed in two-dimensional quantum mechanics [14]. From Eq. (15), it follows that the 
canonical momentum corresponding to the angular variable is 

fij 

Pe = -j- = (r 2 + 8). (16) 
36 

Thus, the classical action / pgdO over an orbit that encloses the origin is shifted by ±2tt8 
due to the existence of the flux line, according to the Aharonov-Bohm effect. However, 
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since the extra term due to the flux line in the Lagrangean (15) is a total time derivative 
(6 being a constant), the classical Lagrangean equations of motion remain unchanged. The 
Hamiltonian corresponding to Eq. (15) is easily found to be 



H = \\ P 2 r + { -^-^] + V HH (r,6), (17) 
2 r z 



where p r = r is the canonical radial momentum. The corresponding Schrodinger equation 
has to be solved numerically again, and its eigenvalue spectrum is altered non trivially. 
We first point out some systematics of this spectrum which will be used in analyzing both 
the quantum density of states and the NNS distributions. The Hamilton operator in the 
presence of the flux line is given by 

Id 2 Id (-i— - S) 2 
H = k-4-,--4~+ { m 2 ' } + Vhh(t, 9). (18) 



2 dr 2 r dr 



This Hamiltonian is diagonalized as before in a large harmonic oscillator basis, but this time 
the basis includes the flux term explicitly. This is possible because the Hamiltonian (18) 
with Vhh replaced by the harmonic oscillator potential |r 2 can be analytically solved with 
eigenvalues 

Ei°\6) = 2n+\l-6\ + l (19) 

and the corresponding eigenvectors |n, |/ — 6\). Again, basis states with A/ = 3 are mixed 
by the nonlinear term of the HH potential. Thus, an eigenfunction has the general form 

v&M) = £ e ^^ie 5 H (20) 

in terms of the radial parts i?^~ 5 '(r) of the harmonic-oscillator basis functions. The bar over 
the first summation sign in Eq. (20) indicates that only the / values (differing by 3 units) 
within one of the sets I, II or III should be summed over. The Hamiltonian (18) may be 
transformed to the original Hamiltonian without flux line by a gauge transformation on the 
wave function: 

M, ^ q,' = e wq = Y^e^^Y^CniRV^r). (21) 
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Since the gauge-transformed Hamiltonian is now independent of 6, the full spectrum must 
remain unchanged under the mapping 8 — > (1 — 6), as this mapping only leads to a reshuffling 
of the / values (recall that < 8 < 1). In particular, the spectrum at 8 = is exactly the 
same as at 8 = 1, since all / values are shifted by one unit and the global properties of the 
system cannot change when all states are included. Obviously, the same symmetry must be 
found in the level density discussed below. This has been checked numerically and found to 
be correct. 



A. Quantum density of states 

In this subsection, the coarse-grained oscillating quantum level density 8g(E) } defined in 
Eq. (13), will be studied for varying flux strength 8. The global analysis immediately shows 
that it is changed drastically in the presence of the flux line. In order to make a detailed 
analysis, we use two windows corresponding to 7 0SC = 0.25 and 1.0, respectively. 

In Figs. 4 and 5, 8g(E) is plotted versus E with a=0.04 and 0.06, respectively, for the 
three values 6=0, 0.25, and 0.5 of the flux strength, together with their Fourier transforms. 
The smooth part g(E) was obtained with 7=1.4 as before, while 7 OSC =0.25 was used in order 
to see finer details of the oscillating level density which will become evident mostly in the 
Fourier spectra. At a=0.04 (Fig. 4) the Fourier spectrum is better resolved, because here we 
could use a larger energy range. In the absence of the flux line (top figures), the three peaks 
corresponding to the averaged periods of orbits A, B, and C can clearly be distinguished; C 
has the lowest period and B the largest. Note that the relative heights of these three peaks 
correspond almost quantitatively to the amplitudes A\ (2) found numerically [cf. Fig. (4) 
of Ref. [7]] and averaged over the energy, if the amplitude of C is multiplied by its time- 
orientation degeneracy factor 2 but those of A and B are not multiplied by their extra 
symmetry degeneracy factor 3. Thus, the Fourier analysis of the quantum level density 
seems to support our choice to give all three orbits the same symmetry degeneracy factor 
which led to Fig. 3. 
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We now consider the effects of the flux line on the beats in Figs. 4 and 5. Notice that only 
one beat minimum persists at all values of 8 that are shown. This minimum is essentially 
due to the interference of orbits A and B. In the special case (5=0.25 it is the only minimum 
present in the energy range considered, and the peak corresponding to orbit C is missing in 
the Fourier spectrum. A simple explanation of this fact will be offered in the next subsection. 
The peak belonging to orbit C appears again at (5=0.5 along with some of the original beat 
structure. 

A very interesting phenomenon appears when 8g(E) is examined with a stronger coarse- 
graining, obtained with 7 OSC =1.0, that corresponds to an emphasis on periods shorter than 
the fundamental period T ~ 1. Since the periods of the primitive orbits A, B, and C are 
all close to one (at least at 8=0) , we should expect all oscillations of the level density to die 
out. This is, indeed, the case as seen on the top row of Fig. 6. However, for < 8 < 1, 
8g(E) shows a remarkably neat oscillatory behaviour with almost no visible beat pattern. 
The corresponding Fourier spectra are also shown in Fig. 6 and exhibit a clear structure 
around the period T w 1/2 (in units of 2tt/uj) in the time domain. The origin of this 
"frequency doubling" in the quantum-mechanical framework is clearly the presence of the 
singular vortex interaction at the origin. When <5=0, the s-state (/ = 0) components of 
the wavefunctions would be non-zero at the origin. But when 8 ^ 0, a genuine zero in the 
wavefunction is introduced at the origin. The frequency doubling may thus be interpreted as 
a diffraction effect of an incoming wave at the location of the flux line: the transmitted wave 
corresponds to the full period, whereas the reflected wave gives rise to a Fourier component 
with a fraction of this period. Indeed, Fig. 6 seems to indicate two peaks around T ~ 0.5. 
In anticipation of the semiclassical interpretation below, we emphasize that this is a pure 
quantum effect, since a classical particle with zero angular momentum cannot penetrate the 
flux line. Also the fact that the quantum-mechanical spectrum is unaltered by the flux line 
when its strength is 8=1 can not be understood classically. 

The scattering of a free charged particle by a singular flux line has been studied in the 
classic paper by Aharonov and Bohm [15]. These authors have shown that the scattering 
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cross section goes like sm 2 (ir6)/ cos 2 (#/2) (lor ^ 7r), where is the scattering angle. It 
vanishes lor both 6=0 and 6=1. Although the particles in the HH potential are not tree 
waves, the leatures discussed above bear some resemblance to the result ol Aharonov and 
Bohm. In particular, the amplitude of the reflected wave (#=0) having a maximum at 
6=0.5 seems to agree qualitatively with the amplitude of the frequency- doubled peak seen 
in the Fourier spectrum of Fig. 6. Unfortunately, the above expression for the scattering 
cross section does not hold for 6=tt } so that no definite conclusions can be drawn for the 
amplitude of the transmitted wave. 

To illustrate this frequency doubling and its genuineness as a quantum phenomenon, 
consider the simple example of particle in a spherical two-dimensional harmonic oscillator 
potential with a flux line at the origin. (This corresponds to choosing a=0 in the HH 
potential.) The exact level density can be computed analytically from the spectrum E^\6) 
given in Eq. (19) and is given by 
1 



g(E,6) 



(huj) 



E \ 1 + 2 J2 cos(2ttM) cos(2irkE/h. 



to 



k=i 



26huj ^ sm(2irk6) sin(7rM) sin(27r kE /2huj) 

k=l k=l 



(22) 



where we have put in the factors %u explicitly. Notice that in the limits 6=0 and 6=1 , 
only the first term remains that gives the exact level density of the unperturbed harmonic 
oscillator: 



E ( 00 1 
9( E ) = 7tT~Y2 1 + 2 E cos(27rLE/M . (23) 
\ hLO ) I k=i ) 



Note that the oscillatory terms are of the form cos(kS/h) with S(E) = 2ttE/uj. The period 
of the fundamental is T = dS/dE = 2tt/uj. In Eq. (22), the corrections due to the flux 
line have an amplitude of order % relative to the leading term. In particular, in the second 
correction term we see that the lowest harmonic has the frequency 2cu, i.e. a period T /2. 
This corresponds to the wave reflected at the singularity, while the transmitted component 
is contained in the first two terms which have the same fundamental period T as the 
unperturbed oscillator. 
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The fact that the new frequency- doubled components appear as % corrections to the 
leading-order terms introduces a caveat for the naive interpretation of these results in terms 
of classical periodic orbits. In order to include diffraction effects, such as those brought 
about by the singular flux line, it is evidently necessary to introduce quantum corrections 
to the standard periodic orbit theory. 

B. Semiclassical analysis and its limitations 

In our semiclassical analysis at 6=0 in Sect. 2 we found a reasonable interpretation of the 
quantum-mechanical level density in terms of the three isolated periodic orbits with shortest 
periods. Now we attempt a semiclassical interpretation of the case with the flux line turned 
on. As emphasized above, we cannot expect to classically simulate the quantum effects of 
diffraction, but we we can still explain some of the observed features in the level density 
and the Fourier spectrum at the leading order. For example, consider the special case with 
6=0.25. As we see in Fig. 4, the contribution from orbit C seems to disappear at this flux 
strength. To understand this, recall that, unlike A or B, the loop orbit C encloses the flux 
line and therefore its action is shifted by ±2tt6 } as already stated above. For 6 = 0.25, 
this shift is ±7r/2 where the plus and minus signs are associated with clockwise and anti- 
clockwise motion of the particle. The total contribution of the orbit C to 6g(E) } assuming 
the amplitude Ac not to depend on the sign of <5, is then found to vanish: 

A c [cos(S c - <j c 7t/2 + tt/2) + cos(S° c - <j c 7t/2 - tt/2)] = 0. (24) 

In the above equation, S c is the action for the orbit C for 6 = 0, and o c is the Maslov 
index whose precise value is irrelevant for this argument. This destructive interference 
between the two time-reversed orbits C by the Aharonov-Bohm phase is, of course, also 
a quantum phenomenon, but it can be incorporated at the leading order of the periodic 
orbit theory. However, at (5=0.25, we still see in Fig. 4 the two periods in the Fourier 
spectrum corresponding to the periods of A and B. Since these orbits do not enclose the 
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flux line, we may assume that they remain unaffected. However, il we now use Eq. (13) and 
omit the contribution from orbit C, the resulting semiclassical Sg(E) reproduces the beat 
minimum seen in Figs. 4 and 5 for 6=0.25 at approximately the right energy, but fails to 
give the correct form of the envelope. This should be clear from the fact that the orbit A 
required for this beat does not exist classically, but is mimicked quantum-mechanically by 
the transmitted wave. The failure becomes even more dramatic at (5=0.5, as shown in Fig. 7, 
if we use Eq. (13) with the unaffected Gutzwiller amplitudes and Maslov indices obtained 
at (5=0, but add the Aharonov-Bohm phase for the orbit C: 



This limitation of the standard semiclassical analysis calls for a systematic investigation 
of quantum corrections to the leading-order terms provided by the periodic classical orbits. 
Such corrections might arise either from extensions of the saddle-point approximation, or 
from corrections going beyond the leading order which would involve new, not necessarily 
closed orbits. 

IV. EFFECT OF THE MAGNETIC FLUX LINE ON THE NNS DISTRIBUTIONS 

While a study of the quantum density of states of a classically chaotic system unravels its 
global behaviour, the local fluctuations are often analyzed by studying the nearest neighbour 
spacing (NNS) distributions in the eigenspectrum [16]. In general the NNS distribution is 
expected to be of Poisson type for integrable systems while it is of Wigner type for classically 
chaotic systems, though exceptions to this rule are known [17]. In this section the NNS 
distributions for the HH system are analyzed with and without the flux line. 

Recall that the HH Hamiltonian is block diagonal and the eigenvalues can be put in 
three distinct and disconnected sets, namely I, II and III (see section I) due to the three-fold 





3*4,4 cos 
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symmetry of the Hamiltonian. The NNS distributions are therefore analyzed in each of these 
classes separately. The symmetry classes are preserved even in the presence of the flux line, 
since the term corresponding to the coupling of the particle to the flux line is separately 
invariant under rotations. The NNS, for a sequence of eigenvalues {£1, E 2} ...} is defined 
through the following unfolding procedure [17] 



where g(E{) is the full density of states of the system evaluated locally and s 8 - is the di- 
mensionless nearest neighbour spacing. The sequence of the s 8 - thus generated are put into 
various bins to generate the NNS distribution denoted by P(s). Typically about 300 or more 
states are included in the analysis to obtain sufficient statistics. This restricts the range of 
the HH parameter to be a < 0.08, since the number of bound states below the barrier energy 
decreases with increasing a. 

We first consider the case 6=0 where there is no flux line. We restrict here the discussion 
to the case a = 0.06, having assured that one finds the expected signatures for the transition 
from regular to chaotic behaviour as a is increased from 0.04 to 0.08. For a = 0.06, we 
have typically about 400 eigenvalues below the barrier energy £'*=46.3 and we scan the 
energy range 16 < E < 43, omitting the low-energy end. The P(s) distributions are 
shown separately for the sets I, II and III in Figs. 8, 9, and 10, respectively. At <5=0, the 
P(s) distributions for Set II (Fig. 9a) and Set III (Fig. 10a) are identical and have more 
resemblance to Wigner-type distributions than the Poisson type, while for Set I (Fig. 8a) 
they are more of Poisson type. This may be understood, as mentioned earlier, from the fact 
that all the states in Set I have the same three-fold symmetry of the HH Hamiltonian since 
/ is always a multiple of 3. This symmetry is lost in the sets II and III which, however, have 
identical eigenvalues and hence identical P(s) distributions. A similar effect is found in the 
case of a particle moving in a rectangular enclosure in the presence of a flux line, where the 
P(s) distribution for states which share the symmetry of the Hamiltonian (in this case the 
parity) is of Poisson type, whereas for other sets it is more like a Wigner distribution [18]. 




(26) 
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In the presence of the flux line, some nontrivial changes are introduced in the NNS 
distribution. We first summarize the changes and then put forward the reasons for the 
same. The effect of the flux line on the NNS distributions is shown in Figs. 8-10 for 8 =0, 
0.231, 0.333, 0.5, 0.769, and 1.0, labeled a,b,...,f within each figure. It is obvious that there 
are some nontrivial features. Even though the P(s) distributions are close to either Poisson 
of Wigner for 8 = and 8 = 1 they do not have these generic feature for arbitrary values of 
8. However, notice that the P(s) distributions for the set I and set II are related (see Figs. 8 
and 9), 



where the subscript denotes the set of eigenvalues. At the outset the mapping 8 to (1 — 8) 
appears peculiar, but actually there is a hint on the analytic structure of the eigenvalues 
even though its full form cannot be obtained analytically. From Eq. (20) it is obvious that 
the eigenvalues in each set in principle depend on all Vs in that given set. Eq. (21) also 
indicates that these Vs are all shifted by 8 in the presence of the flux line. However we 
still do not know how the eigenvalues depend on / and 8. The above relations between 
P(s) distributions in each set suggest that the eigenvalues can only be functions of the form 
f(\k — 8\) with the dependence of /; and 8 restricted to be of the form, — 8\. To see how 
this works, consider the mapping 8 — > (1 — 8), then 




(27) 



and also (see Fig. 10) 






(29) 



That is a shift in <5 — > (1 — <5) is equivalent to /; 



(li — 1) and therefore 



set I: {...,—6,-3,0,3,6,...} 



<-> set II: {...,7,4, 1,-2,-5,...}, 



set III: {...,-7,-4,-1,2,5,...} 



<-> set III: {...,8,5,2,-1,-4,...}, 



(30) 
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which explains the observed results for 6 — > (1 — 6). In the absence of the cubic term in 
VhHi it is known analytically that the eigenvalues depend on / and 6 through |/ — 6\. It is 
therefore interesting that this form carries over to the full HH potential with the flux line. 

It may also be noticed that the envelope of the P(s) distribution changes from Poisson 
to almost Wigner shape in set I as a function of <5, while the opposite is the case for set II. 
This quantum effect may be understood by noting that as 6 is increased, set II goes over to 
set I, thereby introducing more symmetry (see the discussion after Eq. 5). Although P(s) 
distributions for a billiard in the presence of a singular flux line have been studied before 
[19], their interpolation property with the fractional strength of the flux line is novel. 

V. SUMMARY 

We have interpreted the coarse-grained quantum level density of the Henon-Heiles po- 
tential in terms of the periodic orbit theory. The pronounced beats can be explained by 
the interference of three types of short periodic orbits A, B, and C which have nearly iden- 
tical actions and periods in the energy region well below the barriers. There is, however, 
a problem of counting the orbits. The conventional counting of different isolated primitive 
orbits, taken together with their numerically determined amplitudes in the Gutzwiller trace 
formula, does not seem to agree with the amplitudes found in the Fourier transform of the 
quantum density, unless an extra weight of 3 is given to the orbit C. Including this extra 
factor, the lowest harmonics of the trace formula give an almost quantitative agreement with 
the coarse-grained quantum density. 

In our future studies, we plan to investigate the counting problem in more detail by 
studying other potentials with isolated periodic orbits. 

The inclusion of a magnetic flux line through the origin causes substantial changes in the 
quantum spectrum, affecting the beats in the level density dramatically. A straight-forward 
application of the trace formula, adding the Aharonov-Bohm phase to the action of the 
loop orbit C that surrounds the flux line, fails to yield a satisfactory agreement with the 
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quantum level density (both with and without the extra factor 3 for orbit C), although the 
disappearance of the Fourier peak corresponding to orbit C at the flux strength 8 = 0.25 
can be understood in terms of the Aharonov-Bohm phase. Additionally, quantum scattering 
of the particle by the singular flux line takes place. This creates a new peak in the Fourier 
spectrum that corresponds to a reflected wave, simulating a new classical orbit of about half 
the period of the other leading orbits, while the Fourier signal of the transmitted wave - 
whose corresponding classical full period should no longer exist - still persists. To emphasize 
this point, the analytical expression of the quantum level density for a harmonic oscillator 
potential with a flux line through the origin is presented. It shows that corrections of order 
h with respect to the leading terms arise in the level density when the flux line is turned on. 

It is likely that the calculation of Maslov indices has to be modified in the presence of the 
vector potential that describes the flux line, in order to take account of the difference between 
the mechanical and the canonical momenta. The determination of higher-order quantum 
corrections to the Gutzwiller formula, in the presence of a flux line, poses a challenging 
problem for future research. 

Although no major surprises occurred in studying the nearest-neighbour-spacing dis- 
tributions of the Henon-Heiles spectrum, they were found to exhibit novel interpolation 
properties as the flux strength is varied from zero to unity. 

We are grateful to Stephen Creagh for his invaluable help in understanding and com- 
puting the Maslov indices. Thanks are also due to J. Lefebvre, A. Magner, H. Nishioka 
and S. Reimann for stimulating discussions. This research was supported by the Natural 
Sciences and Engineering Research Council of Canada. M.B. and M.V.N.M. acknowledge 
the hospitality of the Department of Physics and Astronomy at McMaster, and R.K.B. that 
of the Institute of Mathematical Sciences at Madras. M.V.N.M. acknowledges the financial 
support from the Hans-Vielberth-Universitatsstiftung Regensburg. 



20 



VI. APPENDIX: MONODROMY MATRIX AND CALCULATION OF THE 

MASLOV INDICES 



For the calculation of the Maslov indices of the three orbits A, B, and C, we have followed 
the methods given by Creagh et al. [8]. In general, the overall index a appearing in the trace 
formula (1) for a given periodic orbit is a sum of two contributions: 



a = fi + v, 



(31) 



whereof the first part fi is the Maslov index occurring in the semiclassical expression of the 
Green's function G(r, r'; E) [1] and counts the number of conjugate points of a given orbit 
at fixed energy. The second contribution v arises when taking the trace of G(r, r'; E) in 
order to arrive at the level density g(E) (or its oscillating part). Whereas both fi and v may 
depend on the starting point along a periodic orbit, their sum has been shown [8] to be a 
topological invariant. 

Two different procedures are used, depending on whether an orbit is stable or unsta- 
ble. We shall sketch these methods very briefly below and refer to Creagh et al. for their 
derivation. But first we recall the definition of the monodromy matrix. 

We start from a periodic orbit given in terms of the coordinates q = (x, y) and momenta 
p = (p X} p y ) as functions of time: q(t) = q(t + T), p(t) = p(t + T). The stability of the 
orbit is given by the propagation of small perturbations <5q(t), 8p{t) away from the exact 
solution qo(^), Po{t): 

q(t) = q (t) + 5q(t); p(t) = p (t) + 6p(t). (32) 

These perturbations satisfy the linearized equations of motion 

q( } = q( } (33) 

V *p(t) J \ -V„(f ) ) { Sp(t) J ' 

with T-periodic coefficients given by the second derivatives of the potential V(q) along the 
orbit: 



d_ 
dt 
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V (*) = \o(t + T) = V tJ (q (t)) 



d 2 V(q) 
dq t dq 3 



(34) 



q=qo(t) 



The time evolution of the solutions of the linear differential equations (33) is given by the 
matrizant Xit): 



6p(t) ) 



X(t) 



«5p(0) ) 



(35) 



with the initial condition X{0) = 1. The value of Xit) after one period T is called the 
monodromy matrix M: 



M = X(T). 



(36) 



According to the Lyapounov theorem, the (4x4) monodromy matrix M has two pairs of 
eigenvalues which are mutually inverse. For a conservative system, two of them are unity, 
corresponding to small perturbations along the orbit. We can therefore, after a transforma- 
tion to the intrinsic coordinate system of the orbit, write M in the form 



M 



M 








\ 



V 



o f 



(37) 



where M is the (2x2) stability matrix. If the orbit is stable, its eigenvalues lie on the unit 
circle with values e ±tx and < x ^ ^ ■ F° r unstable orbits, the eigenvalues are real of the 
form ±e ±x , x > being the Lyapounov exponent. In our numerical calculations, we have 
obtained M using the monodromy matrix method of Baranger et al. [22]. 

a) Unstable orbits: 

For an unstable orbit it is possible to obtain the total Maslov index a at once from 
the full monodromy matrix M. For that purpose one chooses, at a suitable starting point 
of the orbit (at time t = 0), the eigenvector e = (e X} e y} e Px} e Py ) that corresponds to the 
unstable eigenvalue of Af, and the flow vector of the Hamiltonian, f = (x } y } p X} p y ). One 
then propagates these two vectors once around the orbit using the matrizant Xit): 
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e(t) = X(t)e , 



f(t) = X(t)f 



(38) 



which yields two time dependent vectors e(t) = e 8 (t) and f(t) = fi(t) (i = 1,2,3,4). Next 
one defines two real (2x2) matrices U (t) and V(t) by 



U(t) 



' 'e 1 (t) Mt)^ 



f2(t) ) 



V(t) 



' e 3 (t) Mt) N 



e 4 (<) U(t) j 



(39) 



The Maslov index a is then found as the winding number of the complex quantity C(t) = 
{det[U(t) — iV(t)]} 2 } following its path in the complex plane during one full period T. 
Fig. 11 shows this path for orbit A at the scaled energy e=0.9 where it is unstable; clearly 
the winding number is a a = 5. 

b) Stable orbits: 

For stable orbits it is not easy to single out a vector that is orthogonal to the orbit, and 
the two contributions to a in Eq. (31) have to be computed separately. 

The calculation of // is obtained by a similar procedure as that outlined above for 
a of an unstable orbit. One takes the flow vector f and combines it with the vec- 
tor e = (0, 0, 1, — x/y), assuming that i/ / at the starting point or, equivalently, 
e = (0, 0, — y/x, 1), if x ^ at the starting point. One propagates these two vectors 
along the orbit using Eq. (38) and defines U according to Eq. (39). The Maslov index // is 
then found as the number of zeros which det£7 assumes during one full period, hereby not 
counting the zero that appears trivially at the starting point. 

The index v is given by the sign of the quantity w defined by 

TrM - 2 



w 



(40) 



where b is the upper right element of the stability matrix: 



M 



1 ^ 

a b 



which is also given by 



(41) 
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b = d ? ±{ j = T l (42) 

Hereby, p±(t) and q±(t) are the momentum and coordinate, respectively, perpendicular to 
the orbit. If w is positive, we have v=0 and if w is negative, then v=\. The calculation of 
b (42) in principle requires a transformation to the intrinsic coordinate system of the orbit 
and is not simple in general. However, for specific cases, such as the stable orbits A and C in 
the HH potential, the sign of b can be found relatively easily numerically, when solving the 
equations of motion on a computer and following the time evolution of a slightly perturbed 
periodic orbit on the screen: one gives it a small perpendicular starting velocity at time 
t and determines the sign of the perpendicular coordinate after one revolution around the 
perturbed orbit, i.e. at time to + T. 

Figure 12 illustrates these methods of determining // and v for orbit A at the scaled 
energy e=0.8 where it is still stable. 

Both above recipes for obtaining // and v can also be used for unstable orbits; we have 
done this as a check of the numerical methods and found the sum // + v always to agree 
with the a found as a winding number as described above. Note that in all the above 
computations of <r, // and z/, one may not start at any turning point of a given periodic orbit 
since singularities or discontinuities can occur at turning points in some of the calculated 
quantities. Thus, for the orbits A we started at the origin, and for the orbits B we started 
at the apex (intersection with a symmetry axis of the HH potential). 

For the special case of a stable isolated periodic orbit without turning points, Gutzwiller 
[1] has shown that, by a somewhat different bookkeeping of the phases, an alternative form 
of the trace formula can be given in which the total Maslov index cancels. The contribution 
of this orbit to Sg(E) is then 

-1 OO rp 

-irE 9 ■ (i- m MM/h), (43) 

Jnr ^ 2sm(fcx/2) 

if the eigenvalues of M are e ±tx with < x ^ ^ ■ It is easy to see that for orbit C with 
the Maslov index crc=3 this leads, indeed, to the same contribution to Sg(E) as that found 
from Eqs. (1) and (2). 

24 



REFERENCES 

[1] M. C. Gutzwiller: Chaos in Classical and Quantum Mechanics (Springer- Verlag, New 
York, 1990); J. Math. Phys. 10, 1004 (1969); ibid., 11, 1791 (1970); ibid., 12, 343 
(1971). 

[2] R. Balian and C. Bloch, Ann. Phys. (NY) 69, 76 (1972). 

[3] V. M. Strutinsky, Yad. Fiz. 3, 614 (1966); Nucl. Phys. A 95, 420 (1967); ibid., A 122, 
1 (1968). 

[4] M. Sieber and F. Steiner, Phys. Lett. A 144, 159 (1990). 
[5] M. Henon and C. Heiles, Astr. J. 69, 73 (1964). 

[6] K. T. R. Davies, T.E. Huston, and M. Baranger, Chaos 2, 215 (1992). 

[7] M. Brack, R. K. Bhaduri, J. Law and M. V. N. Murthy, Phys. Rev. Lett. 70, 568 (1993). 

[8] S. C. Creagh, J. M. Robbins and R. G. Littlejohn, Phys. Rev. A42, 1907 (1990). 

[9] M. Brack and H. C. Pauli, Nucl. Phys. A 207, 401 (1973). 

[10] R. K. Bhaduri and C. K. Ross, Phys. Rev. Lett. 27, 606 (1971); B. K. Jennings, Ann. 
Phys. (NY) 84, 1 (1974). 

[11] Note that the fits given in Ref. [7], although numerically valid in the energy regions 
given, did not respect this harmonic-oscillator limit. 

[12] Ch. Maier, Thesis, University of Regensburg (1993, unpublished). 

[13] R. C. Churchill, G. Pecelli, and D. L. Rod, in Stochastic Behaviour in Classical and 
Quantum Hamiltonian Systems, Eds. G. Casati and J. Ford (Springer- Verlag, N.Y., 
1979) p. 76. 

[14] R. Jackiw, Ann. Phys. 201, 83 (1990); A. S. Goldhaber and R. Mackenzie, Phys. Lett. 
214B, 471 (1990). 

25 



[15] Aharonov and D. Bohm, Phys. Rev. 115, 485 (1959). 

[16] S. W. McDonald and A. N. Kaufman, Phys. Rev. Lett. 42, 1189 (1979). 

[17] Hua Wu, M. Vallieres, D. H. Feng, and D. W. L. Sprung, Phys. Rev. A42, 1027 (1990). 

[18] G. Date, S. R. Jain and M. V. N. Murthy, Preprint IMSc/93-47. 

[19] M. V. Berry and M. Rubnik, J. Phys. A: Math. Gen. 19, 649 (1986); M. Rubnik and 
M. V. Berry, J. Phys. A: Math. Gen. 19, 669 (1986). 

[20] J. M. Leinaas and J. Myrheim, Nuovo Cim. 370, 1 (1977). 

[21] S. Reimann, M. Brack, A. Magner and M. V. N. Murthy, to be published. 

[22] M. Baranger, K. T. R. Davies, and J. H. Mahoney, Ann. Phys. (N.Y.) 186, 95 (1988). 



26 



FIGURES 

FIG. 1. Equipotential lines of the Henon-Heires potential in the plane £, r] of scaled coordinates. 
The straight lines correspond to the threshold energy E = E* . The three classical periodic orbits 
of shortest lengths A ("linear"), B ("smiley") and C ("loop") for E = E* (i.e., e = 1) are displayed 
by the heavy lines. 

FIG. 2. Oscillating part of the level density of the Henon-Heiles Hamiltonian, for the values 
a=0.04 (left), 0.06 (middle) and 0.08 (right). Upper row: Coarse-grained quantum-mechanical 
results defined by Eq. (9). Lower row: Result of the classical periodic orbit calculation, Eq. (13). 

FIG. 3. Same as Fig. 2, but in the lower row, the contribution of orbit C to Eq. (13) has been 
multiplied by 3. 

FIG. 4. Left column: Quantum-mechanical level density Sg(E) in the presence of a flux line 
with strength S=0.0 (top), 0.25 (middle), and 0.5 (bottom) for the HH potential with a=0.04. The 
coarse-graining was done with a smoothing width j osc = 0.25. Right column: Fourier Transform 
of the Sg(E) shown in the left column. The s-axis gives the time period T in units of 2tt/uj. 

FIG. 5. Same as Fig. 4, for a = 0.06. 

FIG. 6. Same as Fig. 4, for a = 0.06 and j osc = 1.0. 

FIG. 7. Oscillating part of the HH level density for a=0.04 without flux line (left) and with flux 
line with strength 8 = 0.5 (right). Top row: Quantum-mechanical result. Middle row: Semiclassical 
result Eq. (25). Bottom row: Same as middle row, but multiplying the contribution of orbit C in 
Eq. (25) by 3. 

FIG. 8. P(s) distributions for the eigenvalue set I for 6 = 0, 0.231, 0.333, 0.5, 0.769, and 1.0 at 
a = 0.06. The dashed line shows a Wigner distribution and the dotted line a Poisson distribution. 
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FIG. 9. Same as Figure 8 for the eigenvalue set IL 

FIG. 10. Same as Figure 8 for the eigenvalue set III. 

FIG. 11. Illustration of the Maslov index calculation for the unstable orbit A at scaled energy 
e=0.9. Plotted is the complex quantity C(t) (see text) in the complex plane over one period. (Start 
and end points are indicated by arrows.) The winding number is 5, giving a a = 5. 

FIG. 12. Illustration of the Maslov index calculation for the stable orbit A at scaled energy 
e=0.8, having a a = 5. Upper part: One period of the slightly perturbed orbit A. (The unperturbed 
orbit A, a straight line at polar angle 30 degrees, is shown by the thin dashed line.) Here we start at 
the origin ^=r]=0 at an angle of 32 degrees, thus with a positive perpendicular velocity p_i_(0). After 
one period T, the perturbed orbit (shown by the heavy line) ends with a positive perpendicular 
coordinate gj_(T), giving b > and thus v=\ (since TtMa < 2; see text). Lower pari: The quantity 
detU (see text) as a function of time / over one full period, exhibiting four zeros (with / > to) and 
thus giving ji=A. 
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